C     ---------------------------------------------------------------
C
C     Postprocessing: calculate various cross sections
C
C     computing total and differential extinction and absorption
C     cross-sections, using the formulas given by Mackowski, Proc. R.
C     Soc. Lond. A 433, 599 (1991) and Xu, Appl. Opt. 36, 9496 (1997)
C
C     ---------------------------------------------------------------
      subroutine cross_section(nL, cext, cabs, csca, cprx, ek, drot,
     $     nmax,uvmax, k,r0,fint, atr1, btr1, rsr, rsi, rsx, px, as,bs
     $     ,indpol)      
      implicit double precision(a-h, o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
      parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
      parameter (nrc=4*np*(np+1)*(np+2)/3+np)
      parameter (nij=nLp*(nLp-1)/2)

      integer u, u0, nmax(nLp), uvmax(nLp), ind(nLp)
      double precision k,r0(6,nLp),drot(nrc,nij),csca,cext,cabs,cprx
      double precision cscai(nLp), cexti(nLp),cabsi(nLp), cprxi(nLp)
      double precision rsr(np,nLp),rsi(np,nLp),rsx(np,nLp),px(np,nLp)
      complex*16 atr1(ni0,nij),btr1(ni0,nij),at(nmp),bt(nmp),ek(np,nij)
      complex*16 as(nLp,nmp),bs(nLp,nmp),as1(nLp,nmp),bs1(nLp,nmp)
      complex*16 A,B,A0,B0, ref(nLp)

      common/fnr/fnr(0:2*(np+2))

C     constants
      pih    = dacos(0.d0)
      twopi  = 4.d0*pih
      pione  = 2.d0*pih

C     初始化需要计算的物理量
      cext = 0.d0
      cabs = 0.d0
      csca = 0.d0
      cprx = 0.d0
      do i = 1,nL
         cexti(i) = 0.d0
         cabsi(i) = 0.d0
         cscai(i) = 0.d0
         cprxi(i) = 0.d0
C     the refraction index of metal sphere
         ref(i) = dcmplx(r0(5,i), r0(6,i))
      enddo

c     the scattering cross section
      call trans(nL, r0, nmax, uvmax, fint, atr1, btr1, ek,
     +              drot, as,bs, as1,bs1, ind)


C
C     1. Scattering Cross Sections
C

c$$$C     Without Pr Cross Section:
c$$$      do i=1,nL
c$$$         do imn=1,uvmax(i)
c$$$            at(imn) = as(i,imn) + as1(i,imn)
c$$$            bt(imn) = bs(i,imn) + bs1(i,imn)
c$$$         enddo
c$$$         do n=1,nmax(i)
c$$$            sc = 0.d0
c$$$            iL = n * n
c$$$            do m = -n,n
c$$$               sc =  sc + dble(dconjg(as(i,iL))*at(iL))
c$$$               sc = (sc + dble(dconjg(bs(i,iL))*bt(iL)))
c$$$               iL = iL + 1
c$$$            enddo
c$$$            cscai(i) = cscai(i) + sc
c$$$            csca = csca + sc
c$$$         enddo
c$$$      enddo

C     Together with Pr Cross Section
      do i=1,nL
         do imn=1,uvmax(i)
            at(imn) = as(i,imn)+as1(i,imn)
            bt(imn) = bs(i,imn)+bs1(i,imn)
         enddo
         do n=1,nmax(i)
            n1 = n+1
            n2 = 2*n
            rn = 1.0d0/dble(n*n1)
            p = fnr(n)*fnr(n+2)/fnr(n2+1)/fnr(n2+3)/dble(n1)
            t = fnr(n-1)*fnr(n+1)/fnr(n2-1)/fnr(n2+1)/dble(n)
            sc=0.d0
            temp=0.d0
            do m=-n,n
               iL = n*(n+1)+m
               sc = sc + dble(dconjg(as(i,iL))*at(iL))
               sc = sc + dble(dconjg(bs(i,iL))*bt(iL))
               rm=dble(m)*rn
               A0=rm*bt(iL)
               B0=rm*at(iL)
               if(n.eq.nmax(i)) goto 51
               u=(n+1)*(n+2)+m
               fnp=fnr(n+m+1)*fnr(n-m+1)*p
               A0=A0+fnp*at(u)
               B0=B0+fnp*bt(u)
 51            if(n.eq.1.or.iabs(m).gt.n-1) goto 52
               u=(n-1)*n+m
               fn=fnr(n+m)*fnr(n-m)*t
               A0 = A0+fn*at(u)
               B0 = B0+fn*bt(u)
 52            temp = temp + dble(dconjg(as(i,iL))*A0)
               temp = temp + dble(dconjg(bs(i,iL))*B0)
            enddo
            cscai(i) = cscai(i) + sc
            csca = csca + sc
            cprxi(i) = cprxi(i) + temp
            cprx = cprx + temp
         enddo
      enddo

C
C     2. Extinction Cross Section
C
      do j=1,nL
         cz=dcos(k*r0(3,j))
         sz=dsin(k*r0(3,j))
         cmz=dcmplx(cz,-sz)
         A=0.d0
         B=0.d0
         do n=1,nmax(j)
            m0 = n*n + n + 1
            u0 = n*n + n - 1
            A = A + fnr(2*n+1)*(as(j,m0) + bs(j,m0))
            B = B + fnr(2*n+1)*(as(j,u0) - bs(j,u0))
         enddo
         if(indpol.lt.1) then
            cexti(j) = cexti(j) + dble((A-B)*cmz)
            cext = cext + dble((A-B)*cmz)
         else
            cexti(j) = cexti(j) - dimag((A+B)*cmz)
            cext = cext - dimag((A+B)*cmz)
         endif
      enddo
C
C     3. Absorbing Cross Section
C
      do j=1,nL
         do n=1,nmax(j)
            A=ref(j)*dcmplx(rsr(n,j), -rsi(n,j))
            temp1=-dimag(A)
            A=px(n,j)*(ref(j)*rsx(n,j)
     +          -dcmplx(rsr(n,j),rsi(n,j)))
            temp=cdabs(A)*cdabs(A)
            if(temp.eq.0.d0) then
               dn=0.d0
            else
               dn=temp1/temp
            endif
            A = dcmplx(r0(5,j), -r0(6,j))*dcmplx(rsr(n,j), -rsi(n,j))
            temp1 = -dimag(A)
            A = px(n,j)*(rsx(n,j) - ref(j)*dcmplx(rsr(n,j),rsi(n,j)))
            temp = cdabs(A)*cdabs(A)
            if(temp.eq.0.d0) then
               cn=0.d0
            else
               cn=temp1/temp
            endif
C
            do m=-n,n
               i = n*n+n+m
               temp1 = dn*cdabs(as(j,i))*cdabs(as(j,i))
     +	             + cn*cdabs(bs(j,i))*cdabs(bs(j,i))
C               if(indpol.lt.1) then
                  cabsi(j) = cabsi(j) + temp1
                  cabs = cabs + temp1
C               else
C                  cabsyi(j) = cabsyi(j) + temp1
C                  cabsy = cabsy + temp1
C               endif
            enddo

         enddo
      enddo

C     Need to amplify with pi and k
      cz = k*k
      csca = 2.d0*twopi*csca/cz
      cext = twopi*cext/cz
      cabs = 2.d0*twopi*cabs/cz
      cprx = 2.d0*twopi*cprx/cz

      return 
      end
